Ejercicio de Simulación

Vamos a hacer un ejercicio de simulación para ver como se identifica la compoente estacional.

library(urca)
library(forecast)
library(tseries)

    ‘tseries’ version: 0.10-47

    ‘tseries’ is a package for time series analysis and
    computational finance.

    See ‘library(help="tseries")’ for details.
library(lmtest)
Loading required package: zoo

Attaching package: ‘zoo’

The following objects are masked from ‘package:base’:

    as.Date, as.Date.numeric
library(uroot)
library(fUnitRoots)
Loading required package: timeDate
Loading required package: timeSeries

Attaching package: ‘timeSeries’

The following object is masked from ‘package:zoo’:

    time<-

Loading required package: fBasics

Attaching package: ‘fUnitRoots’

The following objects are masked from ‘package:urca’:

    punitroot, qunitroot, unitrootTable
library(sarima)
Loading required package: FitAR
Loading required package: lattice
Loading required package: leaps
Loading required package: ltsa
Loading required package: bestglm

Attaching package: ‘FitAR’

The following object is masked from ‘package:forecast’:

    BoxCox

Loading required package: stats4
require("PolynomF")
Loading required package: PolynomF
###Simulación de un proceso con raíz unitaria estacional
#x11()
x <- ts(sarima::sim_sarima(n=144, model = list(iorder=0, siorder=1, nseasons=12, sigma2 = 1),n.start=24),frequency = 12)
plot(x)

acf(x,lag.max = 36)

monthplot(x)

nsdiffs(x)####Decreta cuantas diferencias estacional a través de la aplicación de 
[1] 1
###Algunas pruebas de raíces unitarias estacionales.


###diferencia estacional
Dx=diff(x,lag=12,differences = 1)###lag:periodo s.
plot(Dx)

acf(Dx,lag.max = 36)

monthplot(Dx)

nsdiffs(Dx)
[1] 0
####Simulación de un SAR
#x11()
x1 <- ts(sim_sarima(n=144, model = list(ar=c(rep(0,11),0.8)),n.start=24),frequency=12)
plot(x1)

acf(x1,lag.max = 36)

monthplot(x1)

nsdiffs(x1)
[1] 1
p <- polynom(c(1,c(rep(0,11),-0.8)))
solve(p)
 [1] -1.0187693+0.0000000i -0.8822801-0.5093846i -0.8822801+0.5093846i
 [4] -0.5093846-0.8822801i -0.5093846+0.8822801i  0.0000000-1.0187693i
 [7]  0.0000000+1.0187693i  0.5093846-0.8822801i  0.5093846+0.8822801i
[10]  0.8822801-0.5093846i  0.8822801+0.5093846i  1.0187693+0.0000000i
abs(solve(p))
 [1] 1.018769 1.018769 1.018769 1.018769 1.018769 1.018769 1.018769
 [8] 1.018769 1.018769 1.018769 1.018769 1.018769
###Note lo cerca que están la raíces de la no estacionariedad del proceso, por eso
####aunque si bien el proceso es estacionario, notamos hay una cercanía a 
####e tener una compoenete estacional.
####El anterior modelo puede escribirse como:
x2 <- ts(sim_sarima(n=144, model=list(sar=0.8, iorder=0, siorder=0, nseasons=12),n.start=24),frequency = 12)
plot(x2)

acf(x2, lag.max=48)

monthplot(x2)

nsdiffs(x2)
[1] 1

Ejermplo Pasajeros

Vamos a ver como se hace el modelamiento completo de la serie de pasajeros.

Iniciaremos con la transformación Box-Cox y las pruebas de raíces Unitarias


    Augmented Dickey-Fuller Test

data:  lAirpassengers
Dickey-Fuller = -1.3232, Lag order = 10, p-value = 0.8582
alternative hypothesis: stationary
p-value greater than printed p-value

Title:
 Augmented Dickey-Fuller Test

Test Results:
  PARAMETER:
    Lag Order: 12
  STATISTIC:
    Dickey-Fuller: 3.7872
  P VALUE:
    0.99 

Description:
 Wed Nov 25 05:33:13 2020 by user: 


############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression none 


Call:
lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.103864 -0.023652 -0.001455  0.022160  0.126649 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
z.lag.1       0.005439   0.001436   3.787 0.000241 ***
z.diff.lag1  -0.198813   0.072080  -2.758 0.006737 ** 
z.diff.lag2  -0.273730   0.073309  -3.734 0.000292 ***
z.diff.lag3  -0.233607   0.072258  -3.233 0.001589 ** 
z.diff.lag4  -0.293133   0.073926  -3.965 0.000126 ***
z.diff.lag5  -0.206562   0.072058  -2.867 0.004915 ** 
z.diff.lag6  -0.266919   0.071493  -3.734 0.000292 ***
z.diff.lag7  -0.234526   0.071847  -3.264 0.001436 ** 
z.diff.lag8  -0.327393   0.073197  -4.473 1.79e-05 ***
z.diff.lag9  -0.198455   0.073623  -2.696 0.008054 ** 
z.diff.lag10 -0.279931   0.072710  -3.850 0.000192 ***
z.diff.lag11 -0.176122   0.073011  -2.412 0.017394 *  
z.diff.lag12  0.627402   0.072683   8.632 3.34e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.04308 on 118 degrees of freedom
Multiple R-squared:  0.8564,    Adjusted R-squared:  0.8406 
F-statistic: 54.14 on 13 and 118 DF,  p-value: < 2.2e-16


Value of test-statistic is: 3.7872 

Critical values for test statistics: 
      1pct  5pct 10pct
tau1 -2.58 -1.95 -1.62


Title:
 Augmented Dickey-Fuller Test

Test Results:
  PARAMETER:
    Lag Order: 12
  STATISTIC:
    Dickey-Fuller: -1.5325
  P VALUE:
    0.7711 

Description:
 Wed Nov 25 05:33:13 2020 by user: 

p-value smaller than printed p-value

Title:
 Augmented Dickey-Fuller Test

Test Results:
  PARAMETER:
    Lag Order: 2
  STATISTIC:
    Dickey-Fuller: -7.6337
  P VALUE:
    0.01 

Description:
 Wed Nov 25 05:33:13 2020 by user: 

Identificación de la componente ARMA estacional y la componente ARMA ordinaria

####################################
There were 18 warnings (use warnings() to see them)
######Diferencia Estacional(continuación AirPassengers)#######
monthplot(dlAirpassengers)
nsdiffs(dlAirpassengers)
[1] 1
nsdiffs(AirPassengers)
[1] 1
DdlAirpassengers=diff(dlAirpassengers,lag=12)###lag=s
#x11()
par(mfrow=c(2,1))

plot(dlAirpassengers)
plot(DdlAirpassengers)

monthplot(DdlAirpassengers)
nsdiffs(DdlAirpassengers)
[1] 0
##Autocorrelogramas
x11()
acf(DdlAirpassengers)
acf(DdlAirpassengers,lag.max = 48, ci.type='ma')
pacf(DdlAirpassengers,lag.max = 48)

Ajuste del Modelo y Análisis de Residuales

##Ajuste del modelo
###Arima Estacional o SARIMA(p=0,d=1,q=1)x(P=0,D=1,Q=1)s=12 con transformación logaritmica

#Modelo MA(1) estacional
modelo = Arima(AirPassengers, c(0, 1, 1),seasonal = list(order = c(0, 1, 1), period = 12),lambda = 0)
coeftest(modelo)

z test of coefficients:

      Estimate Std. Error z value  Pr(>|z|)    
ma1  -0.401828   0.089644 -4.4825 7.378e-06 ***
sma1 -0.556945   0.073100 -7.6190 2.557e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
modeloalter= Arima(AirPassengers, c(1, 1, 0),seasonal = list(order = c(1, 1, 0), period = 12),lambda = 0)

## Análisis de residuales
#x11()
residuales <- modelo$residuals
plot(residuales)

acf(residuales)

pacf(residuales)
######Análisis de Outliers
#Test de normalidad
jarque.bera.test(residuales)

    Jarque Bera Test

data:  residuales
X-squared = 5.2265, df = 2, p-value = 0.0733
#Test de autocorrelaci?n
Box.test(residuales, lag = (length(residuales)/4), type = "Ljung-Box", fitdf = 2)

    Box-Ljung test

data:  residuales
X-squared = 37.874, df = 34, p-value = 0.2969
###Estad?ticas CUSUM
res=residuales
cum=cumsum(res)/sd(res)
N=length(res)
cumq=cumsum(res^2)/sum(res^2)
Af=0.948 ###Cuantil del 95% para la estad?stica cusum
co=0.14422####Valor del cuantil aproximado para cusumsq para n/2
LS=Af*sqrt(N)+2*Af*c(1:length(res))/sqrt(N)
LI=-LS
LQS=co+(1:length(res))/N
LQI=-co+(1:length(res))/N
par(mfrow=c(2,1))

plot(cum,type="l",ylim=c(min(LI),max(LS)),xlab="t",ylab="",main="CUSUM")
lines(LS,type="S",col="red")
lines(LI,type="S",col="red")
#CUSUM Square
plot(cumq,type="l",xlab="t",ylab="",main="CUSUMSQ")                      
lines(LQS,type="S",col="red")                                                                           
lines(LQI,type="S",col="red")

Pronóstico

#x11()
Pronosticos=forecast(modelo,h=12,level=0.95)
plot(Pronosticos)

predic<-predict(modelo,n.ahead=12)
plot(predic$pred)



#####Comparación de pronósticos####
library(fpp)
train <- window(AirPassengers,start=c(1949,01),end=c(1959,12))
test <- window(AirPassengers,start=c(1960,01),end=c(1960,12))
fitmodelo <- Arima(AirPassengers, c(0, 1, 1),seasonal = list(order = c(0, 1, 1), period = 12),lambda = 0)
refit <- Arima(AirPassengers, model=fitmodelo)
fc <- window(fitted(refit), start=c(1960,1))


h <- 1
train <- window(AirPassengers,start=c(1949,01),end=c(1959,12))
test <- window(AirPassengers,start=c(1960,01),end=c(1960,12))
n <- length(test) - h + 1
fitmodelo <- Arima(AirPassengers, c(0, 1, 1),seasonal = list(order = c(0, 1, 1), period = 12),lambda = 0)
fc <- ts(numeric(n), start=c(1960,01), freq=12)
for(i in 1:n)
{  
  x <- window(AirPassengers, end=c(1959, 12+(i-1)))
  refit <- Arima(x, model=fitmodelo)
  fc[i] <- forecast(refit, h=h)$mean[h]
}
dife=(test-fc)^2
ecm=(1/(length(test)))*sum(dife)
ecm
[1] 342.0955
LS0tCnRpdGxlOiAiU0FSSU1BIgpvdXRwdXQ6IAogIGdpdGh1Yl9kb2N1bWVudDogZGVmYXVsdAogIGh0bWxfbm90ZWJvb2s6IGRlZmF1bHQKLS0tCgpgYGB7ciBzZXR1cCwgaW5jbHVkZT1GQUxTRX0Ka25pdHI6Om9wdHNfY2h1bmskc2V0KGVjaG8gPSBUUlVFKQpgYGAKCiMjIEVqZXJjaWNpbyBkZSBTaW11bGFjacOzbgoKVmFtb3MgYSBoYWNlciB1biBlamVyY2ljaW8gZGUgc2ltdWxhY2nDs24gcGFyYSB2ZXIgY29tbyBzZSBpZGVudGlmaWNhIGxhIGNvbXBvZW50ZSBlc3RhY2lvbmFsLgoKYGBge3IgU2ltdWxhY2lvbn0KbGlicmFyeSh1cmNhKQpsaWJyYXJ5KGZvcmVjYXN0KQpsaWJyYXJ5KHRzZXJpZXMpCmxpYnJhcnkobG10ZXN0KQpsaWJyYXJ5KHVyb290KQpsaWJyYXJ5KGZVbml0Um9vdHMpCmxpYnJhcnkoc2FyaW1hKQpyZXF1aXJlKCJQb2x5bm9tRiIpCgojIyNTaW11bGFjacOzbiBkZSB1biBwcm9jZXNvIGNvbiByYcOteiB1bml0YXJpYSBlc3RhY2lvbmFsCiN4MTEoKQp4IDwtIHRzKHNhcmltYTo6c2ltX3NhcmltYShuPTE0NCwgbW9kZWwgPSBsaXN0KGlvcmRlcj0wLCBzaW9yZGVyPTEsIG5zZWFzb25zPTEyLCBzaWdtYTIgPSAxKSxuLnN0YXJ0PTI0KSxmcmVxdWVuY3kgPSAxMikKcGxvdCh4KQphY2YoeCxsYWcubWF4ID0gMzYpCm1vbnRocGxvdCh4KQpuc2RpZmZzKHgpIyMjI0RlY3JldGEgY3VhbnRhcyBkaWZlcmVuY2lhcyBlc3RhY2lvbmFsIGEgdHJhdsOpcyBkZSBsYSBhcGxpY2FjacOzbiBkZSAKIyMjQWxndW5hcyBwcnVlYmFzIGRlIHJhw61jZXMgdW5pdGFyaWFzIGVzdGFjaW9uYWxlcy4KCgojIyNkaWZlcmVuY2lhIGVzdGFjaW9uYWwKRHg9ZGlmZih4LGxhZz0xMixkaWZmZXJlbmNlcyA9IDEpIyMjbGFnOnBlcmlvZG8gcy4KcGxvdChEeCkKYWNmKER4LGxhZy5tYXggPSAzNikKbW9udGhwbG90KER4KQpuc2RpZmZzKER4KQojIyMjU2ltdWxhY2nDs24gZGUgdW4gU0FSCiN4MTEoKQp4MSA8LSB0cyhzaW1fc2FyaW1hKG49MTQ0LCBtb2RlbCA9IGxpc3QoYXI9YyhyZXAoMCwxMSksMC44KSksbi5zdGFydD0yNCksZnJlcXVlbmN5PTEyKQpwbG90KHgxKQphY2YoeDEsbGFnLm1heCA9IDM2KQptb250aHBsb3QoeDEpCm5zZGlmZnMoeDEpCnAgPC0gcG9seW5vbShjKDEsYyhyZXAoMCwxMSksLTAuOCkpKQpzb2x2ZShwKQphYnMoc29sdmUocCkpCiMjI05vdGUgbG8gY2VyY2EgcXVlIGVzdMOhbiBsYSByYcOtY2VzIGRlIGxhIG5vIGVzdGFjaW9uYXJpZWRhZCBkZWwgcHJvY2VzbywgcG9yIGVzbwojIyMjYXVucXVlIHNpIGJpZW4gZWwgcHJvY2VzbyBlcyBlc3RhY2lvbmFyaW8sIG5vdGFtb3MgaGF5IHVuYSBjZXJjYW7DrWEgYSAKIyMjI2UgdGVuZXIgdW5hIGNvbXBvZW5ldGUgZXN0YWNpb25hbC4KIyMjI0VsIGFudGVyaW9yIG1vZGVsbyBwdWVkZSBlc2NyaWJpcnNlIGNvbW86CngyIDwtIHRzKHNpbV9zYXJpbWEobj0xNDQsIG1vZGVsPWxpc3Qoc2FyPTAuOCwgaW9yZGVyPTAsIHNpb3JkZXI9MCwgbnNlYXNvbnM9MTIpLG4uc3RhcnQ9MjQpLGZyZXF1ZW5jeSA9IDEyKQpwbG90KHgyKQphY2YoeDIsIGxhZy5tYXg9NDgpCm1vbnRocGxvdCh4MikKbnNkaWZmcyh4MikKYGBgCgojIyBFamVybXBsbyBQYXNhamVyb3MKClZhbW9zIGEgdmVyIGNvbW8gc2UgaGFjZSBlbCBtb2RlbGFtaWVudG8gY29tcGxldG8gZGUgbGEgc2VyaWUgZGUgcGFzYWplcm9zLgoKSW5pY2lhcmVtb3MgY29uIGxhIHRyYW5zZm9ybWFjacOzbiBCb3gtQ294IHkgbGFzIHBydWViYXMgZGUgcmHDrWNlcyBVbml0YXJpYXMKCmBgYHtyIHByZXNzdXJlLCBlY2hvPUZBTFNFfQojIyMjIyNBanVzdGUgU2VyaWUgZGF0b3MgQWlyUGFzc2VuZ2VycyBwb3IgbWVkaW8gZGVsCiNtb2RlbG8gQVJJTUEgZXN0YWNpb25hbCwgYXPDrSBjb21vIHN1IGNvcnJlc3BvbmRpZW50ZQojYW7DoWxpc2lzIGRlIHJlc2lkdWFsZXMgeSBwcm9uw7NzdGljb3MKCmRhdGEoQWlyUGFzc2VuZ2VycykKcGxvdChBaXJQYXNzZW5nZXJzKQojIyMjIyMjIyMjRGVzcHXDqXMgZGUgaGFiZXIgYXBsaWNhZG8gbGEgZGlmZXJlbmNpYSBvcmRpbmFyaWEgeSBlc3RhY2lvbmFsCiMjIyMjI1Byb2NlZGVtb3MgYSB0cmF0YXIgZGUgaWRlbnRpZmljYXIgbGEgZXN0cnVjdHVyYSBkZSBhdXRvY29ycmVsYWNpw7NuCiMjI2EgY29ydG8gcGxhem8oQVJNQSkgeSBlc3RhY2lvbmFsIFNBUk1BLgojIyNQYXJhIGVzbywgZXMgbmVjZXNhcmlvIGhhYmVyIGNvbnZlcnRpZG8gbGEgc2VyaWUgYSBlc3RhY2lvbmFyaWEjIyMKbEFpcnBhc3NlbmdlcnM9bG9nKEFpclBhc3NlbmdlcnMpCnBsb3QobEFpcnBhc3NlbmdlcnMpCiMjIyMjUHJ1ZWJhIGRlIERpY2tleSBGdWxsZXIjIyMjIyMKdHNlcmllczo6YWRmLnRlc3QobEFpcnBhc3NlbmdlcnMsaz0xMCkKZlVuaXRSb290czo6YWRmVGVzdChsQWlycGFzc2VuZ2VycyxsYWdzID0gMTIsdHlwZT0nbmMnKSAgICMjI0hheSBsYSBwcmVzZW5jaWEgZGUgUmHDrXogVW5pdGFyaWEKc3VtbWFyeSh1cmNhOjp1ci5kZihsQWlycGFzc2VuZ2VycywgbGFncyA9IDEyKSkKYWRmVGVzdChsQWlycGFzc2VuZ2VycyxsYWdzPTEyLHR5cGU9J2N0JykgICMjIyNQdWVkZSB0YW1iacOpbiBpbmRpY2FyCiMjIyNMYSBwcmVzZW5jaWEgZGUgdW5hIHRlbmRlbmNpYSBkZXRlcm1pbsOtc3RpY2EKCiMjIyNEaWZlcmVuY2lhIE9yZGluYXJpYSMjIyMjIyMjIyMjIwoKZGxBaXJwYXNzZW5nZXJzPWRpZmYobEFpcnBhc3NlbmdlcnMsbGFnPTEpCiN4MTEoKQpwYXIobWZyb3c9YygyLDEpKQpwbG90KGxBaXJwYXNzZW5nZXJzKQpwbG90KGRsQWlycGFzc2VuZ2VycykKCmFkZlRlc3QoZGxBaXJwYXNzZW5nZXJzLGxhZ3MgPSAyLHR5cGU9J25jJykgICMjI05vIHNlIGRlYmUgZGlmZXJlbmNpYXIgbcOhcyMjIwoKYGBgCgojIyBJZGVudGlmaWNhY2nDs24gZGUgbGEgY29tcG9uZW50ZSBBUk1BIGVzdGFjaW9uYWwgeSBsYSBjb21wb25lbnRlIEFSTUEgb3JkaW5hcmlhCgpgYGB7ciBDb21wb25lbnRlIEVzdGFjaW9uYWwgeSBPcmRpbmFyaWF9CiMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIwojIyMjIyNEaWZlcmVuY2lhIEVzdGFjaW9uYWwoY29udGludWFjacOzbiBBaXJQYXNzZW5nZXJzKSMjIyMjIyMKbW9udGhwbG90KGRsQWlycGFzc2VuZ2VycykKbnNkaWZmcyhkbEFpcnBhc3NlbmdlcnMpCm5zZGlmZnMoQWlyUGFzc2VuZ2VycykKCkRkbEFpcnBhc3NlbmdlcnM9ZGlmZihkbEFpcnBhc3NlbmdlcnMsbGFnPTEyKSMjI2xhZz1zCiN4MTEoKQpwYXIobWZyb3c9YygyLDEpKQpwbG90KGRsQWlycGFzc2VuZ2VycykKcGxvdChEZGxBaXJwYXNzZW5nZXJzKQptb250aHBsb3QoRGRsQWlycGFzc2VuZ2VycykKbnNkaWZmcyhEZGxBaXJwYXNzZW5nZXJzKQoKCiMjQXV0b2NvcnJlbG9ncmFtYXMKI3gxMSgpCmFjZihEZGxBaXJwYXNzZW5nZXJzKQphY2YoRGRsQWlycGFzc2VuZ2VycyxsYWcubWF4ID0gNDgsIGNpLnR5cGU9J21hJykKcGFjZihEZGxBaXJwYXNzZW5nZXJzLGxhZy5tYXggPSA0OCkKYGBgCgoKIyMgQWp1c3RlIGRlbCBNb2RlbG8geSBBbsOhbGlzaXMgZGUgUmVzaWR1YWxlcwpgYGB7ciBBanVzdGVzIHkgUmVzaWR1YWxlc30KIyNBanVzdGUgZGVsIG1vZGVsbwojIyNBcmltYSBFc3RhY2lvbmFsIG8gU0FSSU1BKHA9MCxkPTEscT0xKXgoUD0wLEQ9MSxRPTEpcz0xMiBjb24gdHJhbnNmb3JtYWNpw7NuIGxvZ2FyaXRtaWNhCgojTW9kZWxvIE1BKDEpIGVzdGFjaW9uYWwKbW9kZWxvID0gQXJpbWEoQWlyUGFzc2VuZ2VycywgYygwLCAxLCAxKSxzZWFzb25hbCA9IGxpc3Qob3JkZXIgPSBjKDAsIDEsIDEpLCBwZXJpb2QgPSAxMiksbGFtYmRhID0gMCkKY29lZnRlc3QobW9kZWxvKQptb2RlbG9hbHRlcj0gQXJpbWEoQWlyUGFzc2VuZ2VycywgYygxLCAxLCAwKSxzZWFzb25hbCA9IGxpc3Qob3JkZXIgPSBjKDEsIDEsIDApLCBwZXJpb2QgPSAxMiksbGFtYmRhID0gMCkKCiMjIEFuw6FsaXNpcyBkZSByZXNpZHVhbGVzCiN4MTEoKQpyZXNpZHVhbGVzIDwtIG1vZGVsbyRyZXNpZHVhbHMKcGxvdChyZXNpZHVhbGVzKQphY2YocmVzaWR1YWxlcykKcGFjZihyZXNpZHVhbGVzKQojIyMjIyNBbsOhbGlzaXMgZGUgT3V0bGllcnMKI1Rlc3QgZGUgbm9ybWFsaWRhZApqYXJxdWUuYmVyYS50ZXN0KHJlc2lkdWFsZXMpCiNUZXN0IGRlIGF1dG9jb3JyZWxhY2k/bgpCb3gudGVzdChyZXNpZHVhbGVzLCBsYWcgPSAobGVuZ3RoKHJlc2lkdWFsZXMpLzQpLCB0eXBlID0gIkxqdW5nLUJveCIsIGZpdGRmID0gMikKCgojIyNFc3RhZD90aWNhcyBDVVNVTQpyZXM9cmVzaWR1YWxlcwpjdW09Y3Vtc3VtKHJlcykvc2QocmVzKQpOPWxlbmd0aChyZXMpCmN1bXE9Y3Vtc3VtKHJlc14yKS9zdW0ocmVzXjIpCkFmPTAuOTQ4ICMjI0N1YW50aWwgZGVsIDk1JSBwYXJhIGxhIGVzdGFkP3N0aWNhIGN1c3VtCmNvPTAuMTQ0MjIjIyMjVmFsb3IgZGVsIGN1YW50aWwgYXByb3hpbWFkbyBwYXJhIGN1c3Vtc3EgcGFyYSBuLzIKTFM9QWYqc3FydChOKSsyKkFmKmMoMTpsZW5ndGgocmVzKSkvc3FydChOKQpMST0tTFMKTFFTPWNvKygxOmxlbmd0aChyZXMpKS9OCkxRST0tY28rKDE6bGVuZ3RoKHJlcykpL04KcGFyKG1mcm93PWMoMiwxKSkKcGxvdChjdW0sdHlwZT0ibCIseWxpbT1jKG1pbihMSSksbWF4KExTKSkseGxhYj0idCIseWxhYj0iIixtYWluPSJDVVNVTSIpCmxpbmVzKExTLHR5cGU9IlMiLGNvbD0icmVkIikKbGluZXMoTEksdHlwZT0iUyIsY29sPSJyZWQiKQojQ1VTVU0gU3F1YXJlCnBsb3QoY3VtcSx0eXBlPSJsIix4bGFiPSJ0Iix5bGFiPSIiLG1haW49IkNVU1VNU1EiKSAgICAgICAgICAgICAgICAgICAgICAKbGluZXMoTFFTLHR5cGU9IlMiLGNvbD0icmVkIikgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAKbGluZXMoTFFJLHR5cGU9IlMiLGNvbD0icmVkIikKCmBgYAoKIyMgUHJvbsOzc3RpY28KCmBgYHtyIFByb25vc3RpY299CiN4MTEoKQpQcm9ub3N0aWNvcz1mb3JlY2FzdChtb2RlbG8saD0xMixsZXZlbD0wLjk1KQpwbG90KFByb25vc3RpY29zKQpwcmVkaWM8LXByZWRpY3QobW9kZWxvLG4uYWhlYWQ9MTIpCnBsb3QocHJlZGljJHByZWQpCgoKIyMjIyNDb21wYXJhY2nDs24gZGUgcHJvbsOzc3RpY29zIyMjIwpsaWJyYXJ5KGZwcCkKdHJhaW4gPC0gd2luZG93KEFpclBhc3NlbmdlcnMsc3RhcnQ9YygxOTQ5LDAxKSxlbmQ9YygxOTU5LDEyKSkKdGVzdCA8LSB3aW5kb3coQWlyUGFzc2VuZ2VycyxzdGFydD1jKDE5NjAsMDEpLGVuZD1jKDE5NjAsMTIpKQpmaXRtb2RlbG8gPC0gQXJpbWEoQWlyUGFzc2VuZ2VycywgYygwLCAxLCAxKSxzZWFzb25hbCA9IGxpc3Qob3JkZXIgPSBjKDAsIDEsIDEpLCBwZXJpb2QgPSAxMiksbGFtYmRhID0gMCkKcmVmaXQgPC0gQXJpbWEoQWlyUGFzc2VuZ2VycywgbW9kZWw9Zml0bW9kZWxvKQpmYyA8LSB3aW5kb3coZml0dGVkKHJlZml0KSwgc3RhcnQ9YygxOTYwLDEpKQoKCmggPC0gMQp0cmFpbiA8LSB3aW5kb3coQWlyUGFzc2VuZ2VycyxzdGFydD1jKDE5NDksMDEpLGVuZD1jKDE5NTksMTIpKQp0ZXN0IDwtIHdpbmRvdyhBaXJQYXNzZW5nZXJzLHN0YXJ0PWMoMTk2MCwwMSksZW5kPWMoMTk2MCwxMikpCm4gPC0gbGVuZ3RoKHRlc3QpIC0gaCArIDEKZml0bW9kZWxvIDwtIEFyaW1hKEFpclBhc3NlbmdlcnMsIGMoMCwgMSwgMSksc2Vhc29uYWwgPSBsaXN0KG9yZGVyID0gYygwLCAxLCAxKSwgcGVyaW9kID0gMTIpLGxhbWJkYSA9IDApCmZjIDwtIHRzKG51bWVyaWMobiksIHN0YXJ0PWMoMTk2MCwwMSksIGZyZXE9MTIpCmZvcihpIGluIDE6bikKeyAgCiAgeCA8LSB3aW5kb3coQWlyUGFzc2VuZ2VycywgZW5kPWMoMTk1OSwgMTIrKGktMSkpKQogIHJlZml0IDwtIEFyaW1hKHgsIG1vZGVsPWZpdG1vZGVsbykKICBmY1tpXSA8LSBmb3JlY2FzdChyZWZpdCwgaD1oKSRtZWFuW2hdCn0KZGlmZT0odGVzdC1mYyleMgplY209KDEvKGxlbmd0aCh0ZXN0KSkpKnN1bShkaWZlKQplY20KYGBgCgoK